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ABSTRACT 

In this paper, we propose a method to study the nature of resonant relaxation in near- 
Keplerian systems. Our technique is based on measuring the fractal dimension of the 
angular momentum trails and we use it to analyze the outcome of iV-body simulations. 
With our method, we can reliably determine the timescale for resonant relaxation, as 
well as the rate of change of angular momentum in this regime. We find that growth of 
angular momentum is more rapid than random walk, but slower than linear growth. We 
also determine the presence of long term correlations, arising from the bounds on angu- 
lar momentum growth. We develop a toy model that reproduces all essential properties 
of angular momentum evolution. 
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1 INTRODUCTION 

Dynamical systems where the gravitational potenti al is dom- 
inated by a single mass are called near-Keplerian dTremainel 
2005), since the orbits of smaller masses in such systems are 
very close to Keplerian conic sections. A characteristic feature 
of these systems is the existence of various distinct timescales 
over which the dynamical variables change. For positions and 
velocities, this is the dynamical time fy,, ~ (a 3 /GM)'/ 2 , where 
M is the mass of the central object and a is a given orbit's 
semimajor axis. The bound orbits (ellipses) precess over the 
precession time f prec ~ [M/Nm^t^m, where m* is the mass of 
smaller objects and N is the number of these objects within a. 
The timescale for the evolution of energy is the relaxation time, 
filx ~ (M 2 /m^N In A)t^ yn . These timescales form an hierarchy 
filx ^ f prec 2> f dyn and lead to thr ee regimes for angular mo- 
mentum evolution dTremaine 1999). 

For short timescales (t < fdyn), the changes in angular 
momentum has no correlation since the torques felt at differ- 
ent parts of the orbit vary rapidly. For intermediate timescales 
( f dyn t < fprec), where we are effectively averaging over or- 
bits, the changes in angular momentum are correlated, since the 
configuration of the orbits change slowly. For long timescales 
(t > fp re c), the correlation is lost again, since the orbits are ran- 
domized as they precess in different directions at different rates. 
The presence of an intermediate regime with enhanced angu- 
lar inomMvtymT_£vo2y^ion_was first recognized in the literature 
bv lRauch & Tremainel (ll99fl). They called this process resonant 
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relaxation, since it is the result of a near resonance between an- 
gular and radial frequencies of the orbit. It plays a central role in 
various scenarios that are proposed to t ake place in the vicinity 
of supermassive black holes (see, e.g.,|Magorrian & Tremaing 
ll999HAlexandeill2008l) . 

There are two major uncertainties regarding resonant re- 
laxation. The first one is the boundaries, especially the upper 
limit, for the intermediate regime. The timescale over which the 
torques are coherent is evidently related to the precession time 
fprec, but the detailed nature of this relationship is unknown. 
Order of magnitude estimates for the timescales are not suffi- 
cient, since in some cases the lifetimes of the stars in the sys- 
tems under question are comparable to the timescales of the dy- 
namical processes ( for an example at the Galactic centre, see 
iMadigan et"ai1l2009f) . 

The nature of evolution of angular momentum over in- 
termediate timescales is also uncertain. In this regime, since 
the torques are correlated, the evolution of angular momentum 
is not going to be diffusive (AL oc A/ 1 / 2 ). iRauch & Tremaing 
1 1996b seem to suggest that the evolution in this regime is bal- 
lis tic (AL oc At , see th eir Eq.6 and Fig. 1) and this is also adopted 
bv lEilon et ai] (l2009): but their simulations contain too few stars 
to draw conclusions in this regard. Indeed, in this work we find 
that the nature of the angular momentum evolution in the inter- 
mediate regime is between diffusive and ballistic. 

We expose this behaviour and determine the timescales by 
carrying out simplified simulations of near-Keplerian systems 
and a fractal analysis of the evolution of angular momentum. In 
section [2] we give a brief review of fractals and present a new 
method for determining fractal dimension, suitable for trails ob- 
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Figure 1. The Koch curve, obtained by repeatedly transforming a line 
segment. The result of the first transformation is shown as a broken line 
in red, to demonstrate the self similar structure more clearly. Each part 
of the Koch curve resting on the straight parts of this broken line is a 
smaller but an otherwise identical copy of itself. 

tained in numerical simulations. In section [3] we describe our 
simulations and the fractal analysis of the angular momentum 
evolution. In section [4] we demonstrate that such an evolution 
can be mimicked by a simple random walk that retains a lim- 
ited term memory. We discuss the limitations and implications 
of our findings in section[5] 



2 FRACTAL DIMENSION 

Fractals iMandelbrotlll982h are geometric objects that exhibit 
a number of properties that make them suitable for modelling 
physical processes and natural structures. The defining property 
of fractals is that their Hausdorff-Besicovitch dimension (here- 
after fractal dimension, D) exceeds their topological dimension. 
A numbe r of method s for c alculating the fractal dimension is 
given by iMandelbro here we give a brief sketch and 

develop a new one that is suitable for our purposes. 

A fractal is a self similar object; that is, part of it exhibits 
similar properties to the whole, sometimes only in a statistical 
sense. If a self similar object is made up of n copies of itself, 
each of which is smaller (in length) by l/m, then the object 
has fractal dimension D = Inn/ In m. It is trivial to see that this 
leads to correct numbers for self similar Euclidean objects. For 
example, any line segment can be thought to be composed of 
n = 3 identical copies of itself, each of which is smaller by 1/3 
(m — 3), giving D — 1 ; or a rectangular prism can be cut into n = 
8 identical copies of itself, each of which is smaller (in length) 
by 1/2 (m = 2), giving D = 3. We can use this technique to 
calculate the dimension of a well known fractal, the Koch curve 
(Fig.QJ. This curve consists of n = 4 identical copies of itself], 
each of which is smaller by 1 /3, leading to a fractal dimension 
£> = ln4/ln3 ~ 1.262. 

Another way to obtain the fractal dimension is to mea- 
sure the length of a curve in increasing detail. Since as we use 
smaller and smaller rulers, we will be resolving more and more 
detail; the me asured length of th e curve L, is a function of the 
ruler length e lMandelbrotl ll967t) : 

L(e)oc E 1 -°. (1) 

As an example, let us assume that we measure the length of the 
Koch curve by a given ruler and obtain L$. When we reduce the 
ruler length by 1/3, we shall be traversing the smaller copies in 
exactly the same manner we traversed the whole curve. Since 
there are four smaller copies, each of which will be measured to 

1 Koch curve can also be seen as consisting of two copies of itself, each 
of which is smaller by 1/ V3. 



have a length 1 /3 times the whole curve, we have L' = 4/3 x L. 
To state in more general terms, if we modify our ruler length by 
e' -j— (l/m)e, the measured length changes by L' (n/m)L. It 
is easy to see that this scaling is satisfied when L <x e 1 ~ ( ln "/ 1,1 m ) . 
This method of determining the (fractal) dimension of a curve 
readily generalizes t o real life curves, w hich are self similar only 
in a statistical sense dMandelbrotll 19671) . 

For our purposes, it is more convenient to interpret a curve 
as the motion trail of an object. As we check the position of 
the object on this trail at decreasing time intervals, we will be 
resolving more details and the total trail length we calculate is 
going to increase. In other words, the trail length is going to be 
a function of the sampling interval %. For example, if we sample 
the motion of an object on the Koch curve 4 times, we will be 
measuring the trail shown as the broken line in FigureQ] leading 
to an increase in length by 4/3, with respect to going from the 
beginning to the end in one step. In more general terms, when 
the sampling interval is modified by %' <— (l/n)%, the measured 
length becomes L' <— (n/m)L. In this approach, the fractal di- 
mension is determined through the relation 

(2) 

This technique of measuring fractal dimension of a trail can be 
easily applied to the results from dynamical simulations. 



3 SIMULATIONS 

3.1 Simulation Method and Parameters 

The system we simulated has three components. At the cen- 
tre lies a stationary supermassive black hole (SMBH) of mass 
M SMBH = 4 x 10 6 Mq. Around the SMBH, we have A^eid = 
1200 field stars each with mass rageid = 2.5Mq, semi-major 
axes distributed in a powerlaw cusp M cusp (r) °c r 3 / 2 from 
a min = 0.0001 to a max = 0.03 pc, with eccentricities between 
e min = and e max = 0.95, distributed following the distribu- 
tion function g(e) oc e. The final component is the massless test 
stars all with semimajor axis a tes t = 0.01 pc, and eccentricities 
e test = 0.2, 0.75, and 0.85 (12 stars each). Other orbital elements 
(inclination, longitude of the ascending node, argument of peri- 
centre and mean anomaly) of all stars are picked at random. 
The exact values chosen do not matter too much, since over the 
course of the simulation the eccentricities are randomized. We 
chose these parameters to have a system that somewhat resem- 
bles the environment of S-stars observed at our Galactic cen- 
tre. There are large uncertainties re garding the star distributions 
at this environment i Merrit fcoid , and references therein), and 
f dyn "C fpiec •C Mx hierarchy may not even exist there. However, 
this condition would hold for a region around a SMBH that de- 
veloped a Bahcall-Wolf cusp teahcall & Wol3l 197d [l977h . so 
we expect our method would be applicable to such systems. 

The details of the code that is used for the simulations is 
going to be explained in detail elsewhere, here we point out only 
the key features. Both the field stars and the test stars feel the 
potential of the SMBH, including the general relativistic (GR) 
correction that leads to the prograde precession of the orbitfl 
The test stars also feel the individual potential of the field stars, 
which leads to retrograde precession of their orbits and changes 

2 We use the treatment of Saha & Tremainej il992l their eq. 30) 
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Figure 2. The evolution of angular momentum (black) and energy (red) 
for a test star with initial eccentricity eo = 0.75. For comparison, a 
Brownian motion (blue) and a fractional Brownian motion (yellow; 
c = 18, n = 4000) data as described in sec.|4]are also plotted. All curves 
are sampled at the same abscissae, and the ordinates are adjusted to have 
unit variance. 
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Figure 3. Measured trail lengths as a function of duration of the 
sampling intervals for angular momentum trails and the twice-broken 
power-law fits made to them. Sets are artificially offset from each other 
to avoid confusion. The gray lines in the background are to lead the eye 
and have slopes — 1/2. 



in angular momentum and energy. The field stars do not feel the 
individual potential of each other, but instead see the potential of 
a smooth cusp that is consistent with their distribution. This ap- 
proximation decreases the time required fo r force computation 
signifi cantly and was already employed by iRauch & Tremaind 
( 1996) in their /V-body simulations. For the interactions between 
the test stars a nd field stars we use a softening kernel (K2) of 
iDehned J200lh . with a softening length 10 4 pc. The units we 
adopted are G = Mfi e i<j = 3000Mq = 1 pc = 1. 

The extended mass distribution in the background cluster 
precesses the orbits in a retrograde fashion, and the GR effects 
lead to prograde precession. By our choice of parameters, these 
two effects cancel each other for a star with a ~ 0.01 pc and 
e ~ 0.6. Mass precession becomes more effective as semi-major 
axes get larger and eccentricities get smaller, while the GR pre- 
cession has the opposite behaviour. 

We carry out the integration of the orbits u sing a high order 
Runge -Kutta-Nystrom method, discovered by IB lanes & Moanl 
teOOfl their SRKNf,). We split the Hamiltoni an into Keplerian 
and perturbation parts ( Kin oshita et al.M 199 lh to increase effi- 
ciency and avoid spurious precession. This scheme advances the 
Keplerian orbital elements correctly, except for the truncation 
and the roundoff errors (the largest accumulated error per star 
amounts to ~ 10~ 5 over the course of the whole simulation). 
In particular, the evolution of the Runge-Lenz vector does not 
exhibit a linear dr ift as in the c ase of potential energy-kinetic 
energy splitting fautetaU 1995b . Our treatment of the GR per- 
turbation leads to a small error in mean motion, but since we 
are interested in changes that take place over many orbits, this 
error is not important. We use shared adaptive time steps, but 
time- s ymmetrize the integration with the method of Irlut et"al] 
Jl995t) . Our simulations last 30 code units which corresponds 
to a few precession times of the slowest precessing test stars. 

We record the energy and angular momentum of each test 
star throughout the simulation. In Figure [2] we show the evo- 
lution of these quantities for a star with initial eccentricity 
eo = 0.75. Even by eye, it is possible to tell that these quan- 
tities show a different behaviour. For comparison, in this figure 
we also plot the curves generated by Brownian and fractional 
Brownian motion, as described in Section[4] 



3.2 Analysis of Simulation Results 

We measure the length of the angular momentum trail of each 
test star by sampling it at intervals ranging from the full length 
of the simulation down to a fraction of the orbital period. The 
dependence of the total measured trail length on the sampling 
interval is shown in Figure[3]for a few stars with different initial 
eccentricities. Starting from long sampling intervals and moving 
towards shorter ones, a number of features can be observed on 
this figure: 

• For long timescales, the curves do not have the slope — 1/2 
(corresponding to dimension D = 2, the value for Gaussian 
random walks ( lMandelbrotlll982l) ), but are somewhat steeper. 
This is to be expected since the energy of an orbit changes 
through relaxation and this process is much slower; hence the 
angular momentum evolution is bounded unlike a true random 
walk, and cannot be described by simple diffusion, even for long 
timescales. 

• There is a marked transition to a more coherent motion as 
indicated by a decrease in the slope. The point of this transition 
can be determined with reasonable accuracy for a given star, but 
it is not common for all stars. 

• Even though the slope decreases, it never becomes zero; 
hence the evolution of angular momentum is never ballistic 
(AL ^ At). 

• This slope can also be determined with reasonable accu- 
racy for a given star, but varies from star to star. 

• Looking at the eccentricity evolution of the stars reveals 
that the transition point is later for the stars that moved into 
e ~ 0.6 region where precession is slow. This is in harmony 
with the expectation that for more slowly precessing stars, the 
coherent torques last longer. 

• The slope increases again for short timescales, but much 
before the period of the stars is reached. This randomization of 
the torques is a result of the stochastic nature of the processes 
that develop the torques and dominates the shorter timescale 
randomization that would result from orbital motion, at least 
down to the timescales we resolve. 

Our results verify the presence of an intermediate regime, 
where the angular momentum evolution is enhanced. The evo- 
lution of angular momentum in this regime is not as rapid as 
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Figure 4. Measured trail lengths as a function of duration of the sam- 
pling intervals for various generated random walk data sets and the 
twice-broken power-law fits made to them. The set in the middle (green) 
has c = 18, n = 4000. The sets above it have varying c = 12, 6, and the 
sets below it have varying n — 2000, 1000,500 (see the text for expla- 
nation of these parameters). All sets have N = 200000 points. To avoid 
confusion, sets are artificially offset from each other; but to make com- 
parison easier, top four and bottom four sets are redrawn with dashed 
black lines at the top and bottom of the figure with identical offsets. The 
gray lines in the background are to lead the eye and have slopes — 1/2. 

ballistic growth AL <x At, but more rapid than diffusive growth 
AL«: A/ 1 / 2 . This manner of evolution can be seen as the general- 
ization of Brownian mot i on c alled fractional Brownian motion. 
Mandelbr ot & van Nessl (l968) descr ibes various proper ties of 
this motion, along with applications. Mand elbrotl jl982h gives 
further generalizations and methods to produce such curves. 
Even though these approaches are mathematically elegant and 
complete, in the next section we propose a simpler model that 
is easier to attach a physical interpretation. 
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Figure 5. Same as Fig.|4]but for bounded random walks. The slopes for 
long sampling intervals are steeper than — 1/2. 

an intermediate regime with lowered slope, whereas for long 
and short timescales the motion is more randomized. The up- 
per bound of the intermediate regime is determined by the pa- 
rameter n: the break occurs when the sampling interval matches 
n x §f . The other parameter c determines the slope in the inter- 
mediate and short timescale regime. Larger c leads to a more 
coherent motion and a slope closer to 0. 

We also generated similar data with bounded random 
walks. For those, we started with a (vector) variable of unit mag- 
nitude |P(f = 0)| = 1, and whenever a step led to |P| > 1.5, we 
took that step in the opposite direction. This limit roughly cor- 
responds to the limit experienced by an orbit starting with initial 
eccentricity eo = 0.75, ./circular M) ~ 1-5. The results from this 
bounded random walks are shown in Figure [5] exhibiting the 
long term correlations similar to angular momentum evolution. 



4 A SIMPLE TOY MODEL FOR ANGULAR 
MOMENTUM EVOLUTION 

One dimensional Brownian motion for a particle's position P(t) 
can be generated as follows. Let the motion over At consist of 
N steps with equal duration 8t = At/N. We choose an initial 
value Prj an d at each step either increase or decrease the value 
of P by &P. We decide which action to take by generating a 
sequence of random numbers Xi uniformly sampled from the 
interval [—0.5,0.5] and choosing a threshold value q = 0. At 
each step, if X, < q we increase the value of P and otherwise 
decrease it. To make the variance of the motion independent of 
the number of steps we choose 8P oc 1 /\/~N. This scheme leads 
to Brownian motion (Ga ussian random w alk) for small 8;, i.e, 
a large number of steps jFalconer]|2003l Chap. 16). It can be 
extended to a vector variable by letting each component perform 
independent Brownian motions. 

We can introduce correlations between the increments of 
the variable P by using a "repository". For this, we keep track 
of the last n values of X, and let our threshold value be propor- 
tional to their average q = c(Xi) n , where c is some constant. 
Here, n determines the length of the correlations and c deter- 
mines their strength. We generated a few sets of data this way 
and measured the lengths of the resulting trails with differing 
sampling durations. 

The curves generated this way (Fig. [4} show very similar 
characteristics to angular momentum evolution. They exhibit 



5 DISCUSSION 

In this work, we studied the evolution of angular momentum 
in near-Keplerian systems by analyzing the outcome of A-body 
simulations. The simulation method we use incorporates cer- 
tain approximations. Our test stars are massless, so they do not 
cause a back-reaction in the surrounding cluster. Furthermore, 
the field stars all have the same mass. A mass spectrum would 
change the granularity of the background potential, affecting 
the applied torques and possibly the coherence timescale. Fi- 
nally, since the field stars do not see the granularity of their own 
potential, their angular momenta do not evolve. This decreases 
the rate of randomization of the background cluster, since vec- 
tor resonant relaxation can change the orientation of the orbits 
on timescales comparable to mass and GR precession for some 
system j3 ( for a comparison of t hese t imescales for the Galactic 
centre, see lKocsis & Tremainej d2010h ). If the torques on a star 
mainly change because of the rearrangement of the background 
cluster, rather than the reorientation of the orbit, this decrease in 
randomization would alter the rate of angular momentum evolu- 
tion. All these approximations limit the domain of applicability 
of our A'-body simulation approach; however, none of them al- 
ter the essential mechanism by which the angular momentum 
evolves. 



3 We thank an anonymous referee for pointing out this shortcoming. 
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We analyzed this evolution by calculating the fractal di- 
mension of the angular momentum trail. With this method it is 
possible to reliably determine the onset of and the rate of evo- 
lution in different regimes. A key result of our analysis is that 
the evolution of angular momentum is neither diffusive nor bal- 
listic ^nan^^rej>ime 1 _asjy^^ studies 
(e.g. lHopman & AlexanderkO ol lEilon et al.l2009h . This seems 
to contradict the re s ults o f the numerical experiments done by 
iRauch & Tremainel dl996l) . The reason for this discrepancy is 
not clear, but we speculate that it arises from the low number of 
stars used in that work. 

We also developed a toy model that reproduces the features 
of angular momentum evolution. This model has adjustable pa- 
rameters with clear physical interpretations. The relation be- 
tween the appropriate values of these parameters and the phys- 
ical variables requires more detailed and extensive analysis, 
which is currently underway. Apart from studying different ini- 
tial conditions, we also plan to analyze the components of the 
torque parallel and perpendicular to the angular momentum 
separately. The nature of these torques can be ver y different 
jRauch & Tremainelll996l : iGiirkan & Hopmanll2007t) . so a sep- 
arate analysis should lead to a better understanding. 

iKoutsoviannisI < |2002b developed fractional Gaussian noise 
generators (FGNGs) similar to our repository model. In that 
work he compares autoregressive moving average (ARMA) 
models to FGNGs, and finds that ARMA models are inferior 
for describing long-term correlations. Th ese models need to 
be m odified to have long term memory dShumwav & Staffer! 
2000), to describe angular momentum evolution in near- 
Keplerian systems. Alternatively, they can be used to describe 
short term correlations for torques, and long term correlations 
can be introduced by taking physical bounds i nto account, as 
is don e here. After this paper was submitted, Madi gan et alJ 
also submitted a paper containing their analysis of this 
problem with ARMA models. They use ARMA(1,1) model, 
which fixes the value of the autocorrelation function for very 
large timescales to zero (see their figure 4) and hence does not 
lead to any correlations beyond a given time. 

The computer programs used to generate and analyze the 
data are available from the author upon request. 
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